Thalamocortical Connections Mature along a Sensorimotor-Association Axis of Cortical Developmental Heterochronicty
Read in Data for Developmental Characterization
Glasser parcel names
#Glasser regions with corresponding labels, tract names, and mesulam assignments
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) #create tract variable with format thalamus_$hemi_$region_autotrack, no dashes -
#Glasser regions assigned to mesulam zones
glasser.assignments <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360-mesulam_economo_yeo-assignments.csv")
glasser.assignments <- merge(glasser.assignments, glasser.regions, by = "label", sort = F)
glasser.assignments <- glasser.assignments %>% select(label, tract, orig_parcelname, mesulam_assignment)S-A axis
S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.assignments, by = c("orig_parcelname"), sort = F)
S.A.axis$SA.axis.bin <- as.numeric(cut2(S.A.axis$SA.axis, g = 5)) #divide the S-A axis into 5 ranked bins, 72 regions per binNeurosynth term maps
#Neurosynth z-score maps generated by nimare. Nimare uses multilevel kernel density analysis- Chi-square analysis + FDR-correction to use the same procedure as Neurosynth
neurosynth.terms <- read.csv("/cbica/projects/thalamocortical_development/Maps/neurosynth/atl-glasser360_neurosynth.csv") %>% select(-regionName) %>% select(-timing) #neurosynth term meta-analytic scores for 123 terms present in both NeuroSynth and the Cognitive Atlas, ordered in lh --> rh glasser order
neurosynth.termlist <- names(neurosynth.terms) #list of terms (cue willy wonka line "dont you want to know our names?" "cant imagine why it wouldnt matter")
neurosynth.terms$label[1:180] <- glasser.regions$label[181:360] #lh first
neurosynth.terms$label[181:360] <- glasser.regions$label[1:180] #rhSpatial permutation (spin) test parcel rotation matrix
perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins
spin.parcels <- glasser.regions #order of complete set of glasser parcels for spinningDataset-specific diffusion measure spreadsheets
FA.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_FA_finalsample_primary.csv") #generated by sample_construction/PNC/tractmeasures_dfs_PNC.R
FA.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_FA_finalsample_primary.csv") #generated by /sample_construction/HCPD/tractmeasures_dfs_HCPD.RDataset-specific tract lists
tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.R
tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.RDevelopmental statistics
setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/") #gam output dir
#list all FA development files in development results directory for accessing. Files were generated by gam_functions/fit_ageGams.R
files <- list.files(getwd(), pattern='FA', ignore.case=T, full.names = F)
filenames <- c()
#generate variable names to assign files data to
for(name in files){
Rname <- gsub('^.{12}|.{4}$', '', name) #remove ".csv" from end of filename and "development_" from start of filename
filenames <- append(filenames, Rname) #save filenames into a character vector
}
#read in files and assign to variables
for(i in 1:length(filenames)){
Rfilename <- sprintf("%s",filenames[i]) #save filename as an individual string
if(grepl("gameffects", Rfilename)) {
x <- read.csv(files[i], header = TRUE)
x <- merge(x, glasser.assignments, by = "tract")
assign(Rfilename, x)
}
if(!grepl("temporal", Rfilename) & !grepl("gameffects", Rfilename)) {
x <- read.csv(files[i], header = TRUE)
x <- merge(x, S.A.axis, by = "tract")
assign(Rfilename, x)
}
if(grepl("sexeffects", Rfilename)) {
x <- read.csv(files[i], header = TRUE)
x <- merge(x, glasser.assignments, by = "tract")
assign(Rfilename, x)
}
}
rm(x)A Spectrum of Thalamocortical Developmental Trajectories
Cortex-wide thalamic connectivity developmental profiles
Number of significant developmental effects
#Function to calculate the number and percent of thalamocortical connections showing a significant developmental change in a given measure
sig.effects <- function(measure, atlas, dataset){
ageeffects.df <- get(sprintf("gameffects_%s_%s_%s", measure, atlas, dataset))
ageeffects.df$significant <- p.adjust(ageeffects.df$Anova.smooth.pvalue, method = c("fdr")) < 0.05 #fdr-corrected significant connections
sigeffect.totaln <- sum(ageeffects.df$significant) #total number of significant connections
sigeffect.percent <- round(sigeffect.totaln/length(ageeffects.df$significant)*100, 2) #percent of significant connections
print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age-related changes in %s", dataset, sigeffect.totaln, sigeffect.percent, measure))
}PNC
## [1] "In pnc, 201 thalamocortical connections (90.13 percent) show significant age-related changes in FA"
HCPD
## [1] "In hcpd, 180 thalamocortical connections (78.26 percent) show significant age-related changes in FA"
Number of significant age by sex interactions
#Function to calculate the number and percent of thalamocortical connections showing a significant age*sex interaction
sig.effects.sex <- function(measure, atlas, dataset){
sexeffects.df <- get(sprintf("sexeffects_%s_%s_%s", measure, atlas, dataset))
sexeffects.df$significant <- p.adjust(sexeffects.df$GAM.int.pvalue, method = c("fdr")) < 0.05 #fdr-corrected
sigeffect.totaln <- sum(sexeffects.df$significant) #total number of significant connections
sigeffect.percent <- round(sigeffect.totaln/length(sexeffects.df$significant)*100, 2) #percent of significant connections
print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age by sex interactions", dataset, sigeffect.totaln, sigeffect.percent))
}PNC
## [1] "In pnc, 0 thalamocortical connections (0 percent) show significant age by sex interactions"
HCPD
## [1] "In hcpd, 4 thalamocortical connections (1.74 percent) show significant age by sex interactions"
Thalamocortical connection GAM smooth functions
PNC
smoothcentered_FA_glasser_pnc.plot <- merge((smoothcentered_FA_glasser_pnc %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_pnc %>% select(tract, GAM.smooth.partialR2)), by = "tract") #zero-centered developmental smooth functions for RH connections + tract-specific partial R2 for plotting
ggplot(smoothcentered_FA_glasser_pnc.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) +
geom_line(linewidth = .3, alpha = .8) +
scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(0, 0.1), oob = squish) +
theme_minimal() +
labs(x="\nAge", y="Connection FA\n") +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
theme(
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_PNC.pdf", device = "pdf", dpi = 500, width = 2.02, height = 1.65)HCPD
smoothcentered_FA_glasser_hcpd.plot <- merge((smoothcentered_FA_glasser_hcpd %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_hcpd %>% select(tract, GAM.smooth.partialR2)), by = "tract")
ggplot(smoothcentered_FA_glasser_hcpd.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) +
geom_line(linewidth = .3, alpha = .8) +
scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.02, 0.075), oob = squish) +
theme_minimal() +
labs(x="\nAge", y="Connection FA\n") +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
scale_y_continuous(limits = c(-0.02, 0.0138)) +
theme(
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_HCPD.pdf", device = "pdf", dpi = 500, width = 2.02 , height = 1.65)Region-specific GAM smooths
#Function for connection-specific GAM spline + derivative plots
connection_trajectory_plot <- function(measure, atlas, dataset, tract_name, display_color, y_ticks){
#Read in participant-level data and model-level fitted values and derivatives
participant.data.df <- get(sprintf("%s.%s.%s", measure, atlas, dataset))
fittedvalues.df <- get(sprintf("fittedvalues_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
derivatives.df <- get(sprintf("derivatives_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
derivatives.df$significant.derivative[derivatives.df$significant.derivative == 0] <- NA
#Connection spline plot with participant-level data
trajectory.plot <- ggplot(data = participant.data.df, aes(x = age, y = get(tract_name))) +
geom_point(color = alpha(display_color, 0.5), size = .1) +
geom_line(data = fittedvalues.df, aes(x = age, y = fitted), color = display_color, linewidth = 1) +
geom_ribbon(data = fittedvalues.df, aes(x = age, y = fitted, ymin = lower, ymax = upper), alpha = .8, linetype = 0, fill = display_color) +
labs(x="\nAge", y=sprintf("%s %s: %s\n", tract_name, measure, toupper(dataset))) +
theme_classic() +
scale_x_continuous(breaks = c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
scale_y_continuous(breaks = y_ticks) +
theme(
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))
#Significant derivatives (rate of age-related change) plot
derivatives.plot <- ggplot(data = derivatives.df) +
geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
scale_fill_gradient(low = alpha(display_color, 0.2), high = display_color, na.value = "white") +
labs(x="\nAge", fill = "Derivative") +
theme_classic() +
scale_x_continuous(breaks=c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
theme(
legend.position = "none",
axis.title.y = element_blank(),
axis.text = element_blank(),
axis.line = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank(),
legend.text = element_text(size = 6, family = "Arial", color = c("black")))
allplots <- list(trajectory.plot, derivatives.plot)
tract.plot <- cowplot::plot_grid(rel_heights = c(16, 3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)
return(tract.plot)
}Primary motor cortex (4)
Left 4: S-A axis rank = 16, mesulam assignment = idiotypic
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.4, 0.45, 0.5, 0.55))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)## smooth.increase.offset
## 1 16.06784
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.35, 0.4, 0.45, 0.5))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)## smooth.increase.offset
## 1 15.45184
Lateral parietal (PF)
Left PF: S-A axis rank = 248, mesulam assignment = heteromodal
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)## smooth.increase.offset
## 1 18.4531
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)## smooth.increase.offset
## 1 16.49456
Lateral prefrontal (8BL)
Left 8BL: S-A axis rank = 342, mesulam assignment = heteromodal
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)## smooth.increase.offset
## 1 21.21106
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)## smooth.increase.offset
## 1 21.91667
Maturational patterns are highly reproducible
Correlation of thalamocortical connection development effect sizes between datasets
Pearson’s R
gameffects.merged <- merge(gameffects_FA_glasser_pnc, gameffects_FA_glasser_hcpd, by = c("tract", "label", "orig_parcelname", "mesulam_assignment"), suffixes = c("_pnc", "_hcpd"))
cor.test(gameffects.merged$GAM.smooth.partialR2_pnc, gameffects.merged$GAM.smooth.partialR2_hcpd)##
## Pearson's product-moment correlation
##
## data: gameffects.merged$GAM.smooth.partialR2_pnc and gameffects.merged$GAM.smooth.partialR2_hcpd
## t = 16.04, df = 219, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6677655 0.7902850
## sample estimates:
## cor
## 0.7349671
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, gameffects.merged, by = c("label", "orig_parcelname", "tract")) #full set of parcel data in rh --> lh order for spinning. spin test null correlations use complete obs only. Each null correlation correspondence statistic is thus computed on a slightly reduced set of data, akin to a jackknife procedure
perm.sphere.p(x = spin.data$GAM.smooth.partialR2_pnc, y = spin.data$GAM.smooth.partialR2_hcpd, perm.id = perm.id.full, corr.type = "pearson") ## [1] 0
Correlation plot
ggplot(gameffects.merged, aes(x = GAM.smooth.partialR2_pnc, y = GAM.smooth.partialR2_hcpd)) +
geom_point(color = c("#FCAB6A"), size = 0.5) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
scale_x_continuous(limits = c(-0.01, 0.17)) +
scale_y_continuous(limits = c(-0.06, 0.15)) +
labs(x="\nPNC", y="HCPD\n") +
theme(
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.99, height = 1.58)Correlation of thalamocortical connection age of maturation between datasets
Pearson’s R
cor.test(gameffects.merged$smooth.increase.offset_pnc, gameffects.merged$smooth.increase.offset_hcpd)##
## Pearson's product-moment correlation
##
## data: gameffects.merged$smooth.increase.offset_pnc and gameffects.merged$smooth.increase.offset_hcpd
## t = 8.6088, df = 153, p-value = 8.379e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4545288 0.6686745
## sample estimates:
## cor
## 0.5712442
Spatial permutation (spin) based p-value
perm.sphere.p(x = spin.data$smooth.increase.offset_pnc, y = spin.data$smooth.increase.offset_hcpd, perm.id = perm.id.full, corr.type = "pearson") ## [1] 2e-04
Correlation plot
ggplot(gameffects.merged, aes(x = smooth.increase.offset_pnc, y = smooth.increase.offset_hcpd)) +
geom_point(color = c("#9c3a80"), size = 0.4) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
scale_x_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0.05, 0)) +
scale_y_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0, 1)) +
labs(x="\nPNC", y="HCPD\n") +
theme(
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) Correlation of developmental magnitude and age of maturation within dataset
PNC
Spearman’s rho
cor.test(gameffects_FA_glasser_pnc$GAM.smooth.partialR2, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_pnc$GAM.smooth.partialR2 and gameffects_FA_glasser_pnc$smooth.increase.offset
## S = 234341, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8242399
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract"))
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") ## [1] 0
Correlation plot
ggplot(gameffects_FA_glasser_pnc, aes(x = GAM.smooth.partialR2, y = smooth.increase.offset)) +
geom_point(color = c("#FCAB6A"), size = 0.5) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
labs(x="\nAge effect (partial R2)", y="Age of maturation (years)\n") +
theme(
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_Ageofmaturation_PNC.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.58)HCPD
cor.test(gameffects_FA_glasser_hcpd$GAM.smooth.partialR2, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_hcpd$GAM.smooth.partialR2 and gameffects_FA_glasser_hcpd$smooth.increase.offset
## S = 137678, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8161002
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract"))
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") ## [1] 0
Early and Late Maturing Thalamocortical Connections
Early and late maturing areas
#Identify parcels that fall within first and fourth quartiles of the age of maturation variable
matstage.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
matstage.df$maturation_bracket <- "mid" #set default bracket to mid
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset < quantile(matstage.df$smooth.increase.offset, 0.25, na.rm = T)] <- "early" #set first quartile bracket to early
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset > quantile(matstage.df$smooth.increase.offset, 0.75, na.rm = T)] <- "late" #set fourth quartile bracket to late
matstage.df <- rbind(matstage.df, data.frame(label = "rh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey
matstage.df <- rbind(matstage.df, data.frame(label = "lh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey
ggseg(.data = matstage.df, atlas = "glasser", mapping = aes(fill = maturation_bracket, colour=I("black"), size=I(.06))) + scale_fill_manual(values = c("#FEC22F", alpha("#323280", 0.97), "white", "white" ), na.value = "gray90") + theme(legend.position = "none")## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi side region label smooth.increase.offset maturation_bracket
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… 16.9 mid
Neurosynth decoding of age of maturation maps
Compute cognitive term-development map spatial correlations
#Function to calculate the correlation between a neurosynth term z-score map and a thalamocortical development map
developmentmap.neurosynth.decoding <- function(df, dev.metric, term){
df.allregions <- left_join(spin.parcels, df, by = c("label", "tract", "orig_parcelname")) #full set of parcel data in rh --> lh order
term.cor <- cor(df.allregions[,term], df.allregions[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between neurosynth term map and development metric
term.results <- data.frame("term" = term, "term.correlation" = term.cor)
return(term.results)
}Developmental timing decoding by Neurosynth cognitive terms
PNC
gameffects_neurosynth_pnc <- merge(gameffects_FA_glasser_pnc, neurosynth.terms, by = "label")
devpattern.neurosynth.pnc <- ldply(neurosynth.termlist, function(t){
developmentmap.neurosynth.decoding(df = gameffects_neurosynth_pnc, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
arrange(desc(term.correlation))devpattern.neurosynth.pnc <- rbind(slice_max(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation map
devpattern.neurosynth.pnc## term term.correlation
## 1 expectancy 0.4099804
## 2 monitoring 0.3968318
## 3 cognitive_control 0.3822346
## 4 reasoning 0.3317204
## 5 response_inhibition 0.3131374
## 6 strategy 0.3119431
## 7 memory_retrieval 0.2883541
## 8 thought 0.2863203
## 9 memory 0.2799527
## 10 recall 0.2692535
## 11 object_recognition -0.2919385
## 12 visual_perception -0.2896572
## 13 morphology -0.2747329
## 14 consolidation -0.2571732
## 15 gaze -0.2507995
## 16 expertise -0.2367694
## 17 coordination -0.2367162
## 18 perception -0.2249922
## 19 facial_expression -0.2183252
## 20 induction -0.2049576
HCPD
gameffects_neurosynth_hcpd <- merge(gameffects_FA_glasser_hcpd, neurosynth.terms, by = "label")
devpattern.neurosynth.hcpd <- ldply(neurosynth.termlist, function(t){
developmentmap.neurosynth.decoding(df = gameffects_neurosynth_hcpd, dev.metric = "smooth.increase.offset", term = t)}) %>%
arrange(desc(term.correlation))devpattern.neurosynth.hcpd <- rbind(slice_max(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10))
devpattern.neurosynth.hcpd## term term.correlation
## 1 expectancy 0.3933764
## 2 sentence_comprehension 0.3887129
## 3 thought 0.3767177
## 4 cognitive_control 0.3477497
## 5 response_inhibition 0.3345571
## 6 monitoring 0.3258089
## 7 recall 0.3033990
## 8 emotion_regulation 0.2986220
## 9 language 0.2927547
## 10 reasoning 0.2664190
## 11 multisensory -0.3268444
## 12 visual_perception -0.3132163
## 13 object_recognition -0.3088589
## 14 induction -0.2857823
## 15 localization -0.2604849
## 16 mental_imagery -0.2576200
## 17 adaptation -0.2537276
## 18 gaze -0.2474086
## 19 fixation -0.2473204
## 20 navigation -0.2404549
Overlapping terms across datasets
merge(devpattern.neurosynth.pnc, devpattern.neurosynth.hcpd, by = "term") %>% arrange(desc(term.correlation.x)) %>% select(term)## term
## 1 expectancy
## 2 monitoring
## 3 cognitive_control
## 4 reasoning
## 5 response_inhibition
## 6 thought
## 7 recall
## 8 induction
## 9 gaze
## 10 visual_perception
## 11 object_recognition
ggplot(devpattern.neurosynth.pnc, aes(x = term.correlation, y = term, color = term.correlation)) +
geom_segment(aes(y = reorder(term, term.correlation), yend = term, x = 0, xend =
term.correlation), color = "#989898", linewidth = 0.2) +
geom_point(size = 2.2, alpha = 0.85) +
theme_light() +
scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.29, 0.35), oob=squish) +
labs(x = "\nCognitive Term-Thalamocortical Development Correlation", y = "NeuroSynth Term\n") +
scale_x_continuous(limits = c(-0.3,0.42), breaks = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4)) +
scale_y_discrete(
labels = c("expectancy"=expression(bold(expectancy)),
"monitoring"=expression(bold(monitoring)),
"cognitive_control"=expression(bold(cognitive_control)),
"reasoning"=expression(bold(reasoning)),
"response_inhibition"=expression(bold(response_inhibition)),
"thought"=expression(bold(thought)),
"recall"=expression(bold(recall)),
"induction"=expression(bold(induction)),
"object_recognition"=expression(bold(object_recognition)),
"visual_perception"=expression(bold(visual_perception)),
"gaze"=expression(bold(gaze)),
"induction"=expression(bold(induction), parse=TRUE))) +
theme(
legend.position = "none",
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_line(color = c("gray90")),
axis.ticks.y = element_blank(),
panel.border = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.ticks = element_line(linewidth = .2))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Neurosynth_Maturationdecoding.pdf", device = "pdf", dpi = 500, width = 3.35, height = 3.07)Term overlap null analysis
#Create 10,000 spins of the HCPD age of maturation map for a neurosynth term overlap null analysis. Each spatially rotated map will be correlated with neurosynth term maps as above and the 10 most positively and negatively correlated terms will be identified. The number of overlapping terms with the empirical PNC terms will then be computed, creating a term overlap null based on spun maps for computing a pspin value.
spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd) #complete data for spinning
x <- spin.data$smooth.increase.offset #specific cortical map to spin
#spin the empirical data
perm.id <- perm.id.full
nroi = dim(perm.id)[1] #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = array(NA,dim=c(nroi,nperm)) #initialize
for (r in 1:nperm) { #spin it
for (i in 1:nroi) {
x.perm[i,r] = x[perm.id[i,r]]
}
}
x.spatialperm <- as.data.frame(x.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(x.perm))))
x.spatialperm$label <- spin.data$label
#compute neurosynth correlations for all spins and identify overlapping terms
compute.termoverlap <- function(perm){
this.perm <- sprintf("perm%s", perm)
nullmap.data <- x.spatialperm %>% select(label, all_of(this.perm))
nullmap.neurosynth <- left_join(nullmap.data, neurosynth.terms, by = "label")
nullmap.neurosynth <- left_join(nullmap.neurosynth, glasser.regions, by = "label")
#null neurosynth correlations for this spin
devpattern.neurosynth.null <- ldply(neurosynth.termlist, function(t){
developmentmap.neurosynth.decoding(df = nullmap.neurosynth, dev.metric = this.perm, term = t)}) %>%
arrange(desc(term.correlation))
#identify top and bottom most correlated terms for this spin
devpattern.neurosynth.null <- rbind(slice_max(devpattern.neurosynth.null, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.null, order_by = term.correlation, n = 10))
#compute overlap
pnc.terms.empirical <- devpattern.neurosynth.pnc %>% select(term)
null.terms <- devpattern.neurosynth.null %>% select(term)
num.overlap <- length(intersect(pnc.terms$term, null.terms$term))
return(num.overlap)
}
term.overlap.null <- lapply(1:nperm, compute.termoverlap)
term.overlap.null <- do.call(rbind, term.overlap.null) %>% as.data.frame %>% set_names("overlapping.terms")
write_rds(term.overlap.null, "/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/termoverlap_null.rds")term.overlap.null <- readRDS("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/termoverlap_null.rds")
sprintf("In %s out of 10,000 spins, the term overlap number was equal to or greater than 11", sum(term.overlap.null$overlapping.terms >= 11))## [1] "In 36 out of 10,000 spins, the term overlap number was equal to or greater than 11"
## [1] 0.0036
Thalamocortical Connectivity Maturation is Organized by the Sensorimotor-Association Axis
gameffects_FA_glasser_pnc <- merge(gameffects_FA_glasser_pnc, (S.A.axis %>% select(tract, SA.axis)), by = "tract")
gameffects_FA_glasser_hcpd <- merge(gameffects_FA_glasser_hcpd, (S.A.axis %>% select(tract, SA.axis)), by = "tract")Age-resolved developmental change plots
PNC
derivatives_FA_glasser_pnc$significant.derivative[derivatives_FA_glasser_pnc$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)
ggplot() +
geom_line(data = derivatives_FA_glasser_pnc, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180) +
scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
scale_y_continuous(expand = c(0.03, 0)) +
theme_classic() +
ylab("Sensorimotor-Association Axis Rank\n") +
xlab("\nAge") +
theme(
legend.position = "none",
axis.line = element_line(size = .2),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.ticks = element_line(linewidth = .2))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Derivativeplot_pnc.pdf", device = "pdf", dpi = 500, width = 2.3, height = 3)HCPD
derivatives_FA_glasser_hcpd$significant.derivative[derivatives_FA_glasser_hcpd$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)
ggplot() +
geom_line(data = derivatives_FA_glasser_hcpd, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180) +
scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
scale_y_continuous(expand = c(0.03, 0)) +
theme_classic() +
ylab("Sensorimotor-Association Axis Rank\n") +
xlab("\nAge") +
theme(
legend.position = "none",
axis.line = element_line(size = .2),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.ticks = element_line(linewidth = .2))Thalamocortical connections mature at progressively older ages across the S-A axis
PNC
Spearman’s rho
cor.test(gameffects_FA_glasser_pnc$SA.axis, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_pnc$SA.axis and gameffects_FA_glasser_pnc$smooth.increase.offset
## S = 683025, p-value = 2.394e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4877183
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract"))
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") #fdr-correct this for axis-based analyses below## [1] 0.00075
Correlation plot
ggplot(gameffects_FA_glasser_pnc, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
geom_point(size = 0.5) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
theme_classic() +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
labs(x="\nS-A axis", y="Age of maturation (years)\n") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_agematuration_pnc.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.62)Brain map
Thalamocortical Age of Maturation
agemat.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
agemat.df <- left_join((spin.data %>% select(label)), agemat.df, by = "label")
agemat.df$cortex <- "cortex"
agemat.df <- rbind(agemat.df, data.frame(label = "rh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey
agemat.df <- rbind(agemat.df, data.frame(label = "lh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey
ggplot() +
geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = smooth.increase.offset, colour=I("black"), size=I(.05))) +
theme_void() +
scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint = 17.5, na.value = "gray93", limits = c(16, 19), oob = squish) +
new_scale_fill() +
geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) +
scale_fill_manual(values = c(alpha("white", 0), "white")) +
theme(legend.position = "none")ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Ageofmaturation_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 5, height = 6.5) Sensorimotor-Association Axis
S.A.axis.plot <- S.A.axis[S.A.axis$tract %in% tractlist.pnc$tract,] %>% select(label, SA.axis)
ggplot() +
geom_brain(data = S.A.axis.plot, atlas = glasser, mapping = aes(fill = SA.axis, colour=I("black"), size=I(.05))) +
theme_void() +
scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint = 180, na.value = "gray93") +
new_scale_fill() +
geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) +
scale_fill_manual(values = c(alpha("white", 0), "white")) +
theme(legend.position = "none")ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 4.7, height = 6.2) HCPD
Spearman’s rho
cor.test(gameffects_FA_glasser_hcpd$SA.axis, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_hcpd$SA.axis and gameffects_FA_glasser_hcpd$smooth.increase.offset
## S = 367590, p-value = 2.934e-12
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.5090025
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract"))
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") #fdr-correct this for axis-based analyses below## [1] 0.00155
Correlation plot
ggplot(gameffects_FA_glasser_hcpd, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
geom_point(size = 0.5) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_color_gradient2(low = "goldenrod1", mid = "#ece4f2", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
theme_classic() +
scale_y_continuous(limits = c(10.7, 22), breaks = c(12, 14, 16, 18, 20, 22)) +
labs(x="\nS-A axis", y="Age of maturation (years)\n") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) Thalamocortical maturation along anatomical axes
#Glasser parcel x,y,z coordinates for defining A-P, M-L, and D-V axes
load(file = "/cbica/projects/thalamocortical_development/Maps/parcellations/surface/hcp_mmp1.0.rda") #RDA file from the brainGraph package https://github.com/cwatson/brainGraph/blob/master/data/hcp_mmp1.0.rda with glasser region x, y, and z MNI coordinates. regions are in lh --> rh order
hcp_mmp1.0$x.mni[1:180] <- hcp_mmp1.0$x.mni[1:180]*-1 #convert R-> L hemisphere coords into medial lateral coords
hcp_mmp1.0$label <- NA
hcp_mmp1.0$label[1:180] <- glasser.regions$label[181:360]
hcp_mmp1.0$label[181:360] <- glasser.regions$label[1:180]
hcp_mmp1.0 <- hcp_mmp1.0 %>% select(-hemi)#C-Mt files for defining the group-level core-matrix gradient
CPt.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractaverage_CPt_primary.csv") %>% select("tract", "label", "orig_parcelname", "mean_CPt")
CPt.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractaverage_CPt_primary.csv") %>% select("tract", "label", "orig_parcelname", "mean_CPt")
CPt.gradient <- merge(CPt.glasser.pnc, CPt.glasser.hcpd, by = c("label", "tract", "orig_parcelname"), all.x = TRUE, all.y = TRUE, suffixes = c("_pnc", "_hcpd"))
CPt.gradient$CPt <- CPt.gradient %>% select(mean_CPt_pnc, mean_CPt_hcpd) %>% rowMeans(na.rm = T)
CPt.gradient <- CPt.gradient %>% select(CPt, label)
hcp_mmp1.0 <- left_join(hcp_mmp1.0, CPt.gradient, by = "label")#Function to compute the correlation between a development map and an anatomical axis (x,y,z) as well as test the significance of the difference between the developmental correlation with the anatomical axis versus the S-A axis
compute_axis_correlation <- function(df.dev, dev.metric, axis){
df.dev.axes <- merge(df.dev, hcp_mmp1.0, by = "label")
df.dev.axes.spin <- left_join(spin.parcels, df.dev.axes, by = c("label", "tract", "orig_parcelname"))
#S-A axis - development correlation (as computed above) for comparison of two overlapping correlations based on dependent groups
SA.axis.cor <- cor(df.dev.axes$SA.axis, df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman"))
#Anatomical axis - development correlation
anatomical.axis.cor <- cor(df.dev.axes[,axis], df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between anatomical axis and development metric
anatomical.axis.pvalue <- perm.sphere.p(x = df.dev.axes.spin[,axis], y = df.dev.axes.spin[,dev.metric], perm.id = perm.id.full, corr.type = "spearman") #spin based p-value for the correlation
#S-A axis - anatomical axis correlation
SA.anatomical.cor <- cor(df.dev.axes$SA.axis, df.dev.axes[,axis], use = "complete.obs", method = c("spearman"))
#Test for significant difference between two correlations based on dependent groups (i.e., correlations with an overlapping input)
cocor.result <- cocor.pvalue <- cocor.dep.groups.overlap(
r.jk = SA.axis.cor,
r.jh = anatomical.axis.cor,
r.kh = SA.anatomical.cor,
n = nrow(df.dev.axes), alternative = "two.sided", test = "hittner2003")
cocor.pvalue <- cocor.result@hittner2003[["p.value"]]
comparison.results <- data.frame(axis = axis, axis.cor = anatomical.axis.cor, axis.cor.pvalue = anatomical.axis.pvalue, SA.cor = SA.axis.cor, cocor.pvalue = cocor.pvalue)
return(comparison.results)
}PNC
Anterior-posterior axis
compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "y.mni")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 y.mni 0.3101209 0.0621 0.4877183 0.0001552829
Medial-lateral axis
compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "x.mni")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 x.mni 0.03024772 0.43765 0.4877183 3.832011e-08
Dorsal-ventral axis
compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "z.mni")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 z.mni 0.1173292 0.26005 0.4877183 0.0002034701
Core-matrix gradient
compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "CPt")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 CPt 0.2897652 0.0287 0.4877183 0.0006272672
FDR-correct axis correlation ps
axis.cor.ps <- c(0.00075, 0.0621, 0.43765, 0.26005, 0.0287) #S-A, A-P, M-L, D-V, C-M
p.adjust(axis.cor.ps, method = c("fdr"))## [1] 0.0037500 0.1035000 0.4376500 0.3250625 0.0717500
FDR-correct correlation comparison ps
cocor.ps <- c(0.0001552829, 3.832011e-08, 0.0002034701, 0.0006272672)
p.adjust(cocor.ps, method = c("fdr"))## [1] 2.712935e-04 1.532804e-07 2.712935e-04 6.272672e-04
#Put above results into df for plotting
axis.results <- data.frame(Axis = factor(), corr = double())
axis.results <- axis.results %>% add_row(Axis = "S-A", corr = 0.4877183)
axis.results <- axis.results %>% add_row(Axis = "A-P", corr = 0.3101209)
axis.results <- axis.results %>% add_row(Axis = "C-M", corr = .2897652)
axis.results <- axis.results %>% add_row(Axis = "D-V", corr = 0.1173292)
axis.results <- axis.results %>% add_row(Axis = "M-L", corr = 0.03024772)
axis.results$Axis <- factor(axis.results$Axis, ordered = TRUE, levels = c("S-A", "A-P", "C-M", "D-V", "M-L"))
ggplot(axis.results, aes(x = Axis, y = corr)) +
geom_segment(aes(x = Axis, xend = Axis, yend =
corr, y = 0), linewidth = 0.2) +
geom_point(aes(fill = Axis), size = 3.6, shape = 22, color = alpha("white", 0)) +
scale_y_continuous(limits = c(-0.01, .5)) +
labs(x="Neuroaxis") +
labs(y="\nAxis Alignment") +
theme_classic() +
scale_fill_manual(values=c("#672975", "#9898BF", "#9898BF", "#9898BF", "#9898BF")) +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.text.x = element_text(vjust = 0.7, angle = 20),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Anatomicalaxes_agematuration_comparison_lollyplot.pdf", device = "pdf", dpi = 500, width = 1.36, height = 1.05)HCPD
Anterior-posterior axis
compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "y.mni")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 y.mni 0.441284 0.02195 0.5090025 0.124339
Medial-lateral axis
compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "x.mni")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 x.mni 0.1124713 0.3695 0.5090025 5.601018e-07
Dorsal-ventral axis
compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "z.mni")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 z.mni 0.0814776 0.38565 0.5090025 1.390815e-05
Core-matrix gradient
compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "CPt")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 CPt 0.3809388 0.0163 0.5090025 0.02169413
Correct axis correlation ps
axis.cor.ps <- c(0.00155, 0.02195, 0.3695, 0.38565, 0.0163) #S-A, A-P, M-L, D-V, C-M
p.adjust(axis.cor.ps, method = c("fdr"))## [1] 0.00775000 0.03658333 0.38565000 0.38565000 0.03658333
Correct correlation comparison ps
cocor.ps <- c(0.124339, 5.601018e-07, 1.390815e-05, 0.02169413)
p.adjust(cocor.ps, method = c("fdr"))## [1] 1.243390e-01 2.240407e-06 2.781630e-05 2.892551e-02
Thalamocortical Connectivity Maturation Mirrors Timescales of Cortical Plasticity
Plasticity marker development data
#Fluctuation amplitude age of decrease onset data from Sydnor et al., 2023, Nat Neurosci https://www.nature.com/articles/s41593-023-01282-y
boldamplitude.agedecrease.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/spatiotemp_dev_plasticity/cortical_maps/AgeofDeclineOnset_FirstNegDeriv.pscalar.nii")
boldamplitude.dev <- data.frame(orig_parcelname = names(boldamplitude.agedecrease.cifti$Parcel), amplitude.agedecline = boldamplitude.agedecrease.cifti$data)#T1/T2 ratio rate of developmental change from Baum et al., 2022, J Neuro https://www.jneurosci.org/content/42/29/5681
myelin.maxdev.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin.roc.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_mean_posterior_annualized_roc_myelin.pscalar.nii") #annualized rate of change
myelin.dev <- data.frame(orig_parcelname = names(myelin.maxdev.cifti$Parcel), myelin.agemaxdev = myelin.maxdev.cifti$data, myelin.annualroc = myelin.roc.cifti$data)#Generate the E:I ratio-age slope map in Glasser HCP-MMP parcels from Zhang, Larsen et al., bioRxiv, https://www.biorxiv.org/content/10.1101/2023.06.22.546023v1.full
library(ciftiTools)
#Calculate E:I age slopes in Yan100 parcels
EI.group.ages <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/age_Yan.csv", header = F) %>% as.data.frame %>% set_names("age") %>% mutate(age = age/12) #average age of 29 age groups
EI.group.ratio <- t(read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_Yan.csv", header = F)) %>% as.data.frame() #average E:I ratio in all of the Yan100 parcels in the 29 age groups
Yan100.EI.ageslopes <- sapply(EI.group.ratio, function(x) lm(x ~ EI.group.ages$age)$coefficients[2]) %>% as.data.frame() %>% set_names("slope") #linear model coefficients for the association between age and E:I ratio in all Yan100 parcels
#Reorder Yan100 parcels to match the released version of the atlas
Yan100.region.reordering <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/Yanatlas_regionmapping.csv") #mapping of regions between version
Yan100.EI.ageslopes$mapping <- Yan100.region.reordering$parcelnumber.released
Yan100.EI.ageslopes <- Yan100.EI.ageslopes %>% arrange(mapping)
#Map parcel-level Yan100 data to fslr 32k vertex-level data
Yan100.densecifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/100Parcels_Yeo2011_17Networks.dscalar.nii")
Yan100.lh.vertices <- Yan100.densecifti$data$cortex_left %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.lh.vertices.slope <- left_join(Yan100.lh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
Yan100.rh.vertices <- Yan100.densecifti$data$cortex_right %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.rh.vertices.slope <- left_join(Yan100.rh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
lh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.leftcortex.csv", header = F)
rh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.rightcortex.csv", header = F)
Yan100.EI.ageslopes.densecifti <- as_cifti(cortexL = Yan100.lh.vertices.slope$slope, cortexR = Yan100.rh.vertices.slope$slope, cortexL_mwall = lh.medialwall$V1, cortexR_mwall = rh.medialwall$V1) #ciftify me captain
write_cifti(Yan100.EI.ageslopes.densecifti, "/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii")
#Parcellate the fslr 32k vertex-level E:I-age slope with the glasser HCP-MMP atlas
command = "-cifti-parcellate /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii /cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser_space-fsLR_den-32k_desc-atlas.dlabel.nii COLUMN /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii -only-numeric"
ciftiTools::run_wb_cmd(command, intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)
detach("package:ciftiTools", unload=TRUE)EI.ageslope.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii")
EIratio.dev <- data.frame(orig_parcelname = names(EI.ageslope.cifti$Parcel), EI.ageslope = EI.ageslope.cifti$data)df.list <- list(boldamplitude.dev, myelin.dev, EIratio.dev) #dfs to merge
plasticity.dev <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list)
plasticity.dev <- left_join(spin.parcels, plasticity.dev, by = "orig_parcelname")Fluctuation amplitude decrease onset map
plasticity.dev$agedecline.protracted <- plasticity.dev$amplitude.agedecline
plasticity.dev$agedecline.protracted[plasticity.dev$agedecline.protracted == "NaN"] <- 23 #lets plot regions with most protracted functional expression of plasticity (amplitude never declines!)
snr.exclude <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/SNRmask_glasser360.csv") %>% filter(SNR.mask == 0) #but not those we excluded for low SNR
plasticity.dev$agedecline.protracted[plasticity.dev$label %in% snr.exclude$label] <- NA
ggplot() +
geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = agedecline.protracted, colour=I("black"), size=I(.03))) +
theme_void() +
scale_fill_gradientn(colors = c(alpha("#323280", 0.2), alpha("#323280", 0.4), alpha("#323280", 0.6), alpha("#323280", 0.9), alpha("#323280", 1)), limits = c(10, 17), oob = squish, na.value = "white") ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry orig_parcelname
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… EMPTY L_10pp_ROI
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## # myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## # agedecline.protracted <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label orig_parcelname tract amplitude.agedecline
## 396 lh_L_10pp L_10pp_ROI thalamus_L_10pp_autotrack 17.18593
## myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396 16.37678 0.008886078 -0.03484702 17.18593
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Boldamplitude_agedecline_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.7, height = 5.4) ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry orig_parcelname
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… EMPTY L_10pp_ROI
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## # myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## # agedecline.protracted <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label orig_parcelname tract amplitude.agedecline
## 396 lh_L_10pp L_10pp_ROI thalamus_L_10pp_autotrack 17.18593
## myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396 16.37678 0.008886078 -0.03484702 17.18593
T1/T2 annualized rate of increase map
ggplot() +
geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = myelin.annualroc, colour=I("black"), size=I(.03))) +
theme_void() +
scale_fill_gradientn(colors = c(alpha("#323280", 1), alpha("#323280", 0.8), alpha("#323280", 0.6), alpha("#323280", 0.4), alpha("#323280", 0.2)), limits = c(0.012, 0.016), oob = squish, na.value = "white") ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry orig_parcelname
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… EMPTY L_10pp_ROI
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## # myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## # agedecline.protracted <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label orig_parcelname tract amplitude.agedecline
## 396 lh_L_10pp L_10pp_ROI thalamus_L_10pp_autotrack 17.18593
## myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396 16.37678 0.008886078 -0.03484702 17.18593
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2ratio_annualizedroc_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.4, height = 5.13) ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry orig_parcelname
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… EMPTY L_10pp_ROI
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## # myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## # agedecline.protracted <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label orig_parcelname tract amplitude.agedecline
## 396 lh_L_10pp L_10pp_ROI thalamus_L_10pp_autotrack 17.18593
## myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396 16.37678 0.008886078 -0.03484702 17.18593
E:I ratio magnitude of developmental decline map
ggplot() +
geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = EI.ageslope, colour=I("black"), size=I(.03))) +
theme_void() +
scale_fill_gradientn(colors = c(alpha("#323280", 0.2), alpha("#323280", 0.4), alpha("#323280", 0.6), alpha("#323280", 0.9), alpha("#323280", 1)), limits = c(-0.04, -0.035), oob = squish, na.value = "white") ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry orig_parcelname
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… EMPTY L_10pp_ROI
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## # myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## # agedecline.protracted <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label orig_parcelname tract amplitude.agedecline
## 396 lh_L_10pp L_10pp_ROI thalamus_L_10pp_autotrack 17.18593
## myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396 16.37678 0.008886078 -0.03484702 17.18593
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIratio_decline_corticalmap.pdf", device = "pdf", dpi = 500, width = 5, height = 5.08) ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry orig_parcelname
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… EMPTY L_10pp_ROI
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## # myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## # agedecline.protracted <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label orig_parcelname tract amplitude.agedecline
## 396 lh_L_10pp L_10pp_ROI thalamus_L_10pp_autotrack 17.18593
## myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396 16.37678 0.008886078 -0.03484702 17.18593
Thalamocortical maturational patterns align to the spatiotemporal unfolding of plasticity readouts
PNC
plasticity.dev.pnc <- left_join(plasticity.dev, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract"))Thalamocortical connectivity maturation - E:I ratio magnitude of developmental decline
Spearman’s rho
cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$EI.ageslope, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$EI.ageslope
## S = 733425, p-value = 2.329e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4499176
Spatial permutation (spin) based p-value
EI.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") Correlation plot
ggplot(plasticity.dev.pnc, aes(x = EI.ageslope, y = smooth.increase.offset)) +
geom_point(size = 0.5, color = "#A2A2C8") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
theme_classic() +
labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) Thalamocortical connectivity maturation - myelin annualized rate of growth
Spearman’s rho
cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$myelin.annualroc, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$myelin.annualroc
## S = 1929655, p-value = 3.141e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4472777
Spatial permutation (spin) based p-value
myelin.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") Correlation plot
ggplot(plasticity.dev.pnc, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
geom_point(size = 0.5, color = "#A2A2C8") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_x_continuous(limits = c(0.005, 0.025)) +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
theme_classic() +
labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) Thalamocortical connectivity maturation - fluctuation amplitude age of decrease onset
Spearman’s rho
cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$amplitude.agedecline, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$amplitude.agedecline
## S = 712074, p-value = 3.088e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3028346
Spatial permutation (spin) based p-value
amplitude.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman")Correlation plot
ggplot(plasticity.dev.pnc, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
geom_point(size = 0.5, color = "#A2A2C8") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_x_continuous(limits = c(10, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
theme_classic() +
labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) FDR-corrected p-values across correlations
## [1] 0.00020 0.00090 0.02965
## [1] 0.00060 0.00135 0.02965
HCPD
plasticity.dev.hcpd <- left_join(plasticity.dev, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract"))Thalamocortical connectivity maturation - E:I ratio magnitude of developmental decline
Spearman’s rho
cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$EI.ageslope, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$EI.ageslope
## S = 409170, p-value = 9.57e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4534636
Spatial permutation (spin) based p-value
EI.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") Correlation plot
ggplot(plasticity.dev.hcpd, aes(x = EI.ageslope, y = smooth.increase.offset, color = EI.ageslope)) +
geom_point(size = 0.5, color = "#A2A2C8") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
theme_classic() +
labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) Thalamocortical connectivity maturation - myelin annualized rate of growth
Spearman’s rho
cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$myelin.annualroc, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$myelin.annualroc
## S = 1071660, p-value = 7.233e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4314374
Spatial permutation (spin) based p-value
myelin.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") Correlation plot
ggplot(plasticity.dev.hcpd, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
geom_point(size = 0.5, color = "#A2A2C8") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_x_continuous(limits = c(0.005, 0.025)) +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
theme_classic() +
labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) Thalamocortical connectivity maturation - fluctuation amplitude age of decrease onset
Spearman’s rho
cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$amplitude.agedecline, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$amplitude.agedecline
## S = 378150, p-value = 7.236e-08
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4136816
Spatial permutation (spin) based p-value
amplitude.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman") Correlation plot
ggplot(plasticity.dev.hcpd, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
geom_point(size = 0.5, color = "#A2A2C8") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_x_continuous(limits = c(10, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
theme_classic() +
labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) FDR-corrected p-values across correlations
## [1] 0.00010 0.00195 0.00890
## [1] 0.000300 0.002925 0.008900